P3768 简单的数学题

i=1nj=1nijgcd(i,j)\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)

k=1ni=1nj=1nij[gcd(i,j)=k]\sum_{k=1}^n\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=k]

k=1nk3i=1nkj=1ij[gcd(i,j)=1]\sum_{k=1}^nk^3\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor}\sum_{j=1}^{}ij[gcd(i,j)=1]

k=1nk3d=1nkμ(d) d2i=1nkdj=1nkdij\sum_{k=1}^nk^3\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor}\mu(d)~d^2\sum_{i=1}^{\lfloor \frac{n}{kd} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{kd} \rfloor}ij

sum2(n)=i=1ni2,sum3(n)=i=1ni3sum^2(n)=\sum_{i=1}^ni^2 , sum^3(n)=\sum_{i=1}^ni^3

利用小学奥数知识有:

sum2(n)=n(n+1)(2n+1)6   ,   sum3(n)=n2(n+1)24sum^2(n)=\frac{n(n+1)(2n+1)}{6} ~~~ , ~~~ sum^3(n)=\frac{n^2(n+1)^2}{4}

上式为:

k=1nk3d=1nkμ(d) d2sum3(nkd)\sum_{k=1}^nk^3\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor}\mu(d)~d^2sum^3(\frac{n}{kd})

T=kdT=kd , 则原式可继续化简为

T=1nkTk3μ(Tk)(Tk)2sum3(nT)\sum_{T=1}^n \sum_{k|T}k^3 \mu(\frac{T}{k})(\frac{T}{k})^2sum^3(\frac{n}{T})

T=1nsum3(nT) T2kTkμ(Tk)\sum_{T=1}^n sum^3(\frac{n}{T})~T^2 \sum_{k|T}k \mu(\frac{T}{k})

后面是一个卷积形式: idμ=φid*\mu=\varphi

T=1nsum3(nT) T2φ(T)\sum_{T=1}^n sum^3(\frac{n}{T})~T^2 \varphi(T)

现在只需求得 φ×id×id\varphi \times id \times id 的前缀和就可以整除分块了。


f(n)=φ(n)×n×nf(n)=\varphi(n) \times n \times n

g(n)=n×ng(n)=n \times n

h(n)=dnf(d)g(nd)h(n)=\sum_{d|n}f(d)* g(\frac{n}{d})

                         =dnφ(d)×d×dndnd~~ ~~~~~~~~~~~~~~~~~~~~~~~ =\sum_{d|n}\varphi(d) \times d \times d * \frac{n}{d} * \frac{n}{d}

      =n2dnφ(d)=n3~~ ~~~~ =n^2\sum_{d|n}\varphi(d)=n^3

g(1)S(n)=i=1nh(i)d=2ng(d)S(nd)g(1)S(n)=\sum_{i=1}^{n}h(i)-\sum_{d=2}^ng(d)S(\lfloor \frac{n}{d} \rfloor)

S(n)=sum3(n)d=2ng(d)S(nd)S(n)=sum^3(n)-\sum_{d=2}^ng(d)S(\lfloor \frac{n}{d} \rfloor)

这道题就完了。

#include <map>
#include <cstdio>
using namespace std;

const int MAXN = 8000000;
long long n;
int p , inv4 , inv6 , k , prime[ MAXN + 5 ] , phi[ MAXN + 5 ] , f[ MAXN + 5 ];
bool vis[ MAXN + 5 ];

void sieve( ) {
	f[ 1 ] = phi[ 1 ] = 1;
	for( int i = 2 ; i <= MAXN ; i ++ ) {
		if( !vis[ i ] ) {
			prime[ ++ k ] = i;
			phi[ i ] = i - 1;
		}
		for( int j = 1 ; j <= k && 1ll * i * prime[ j ] <= MAXN ; j ++ ) {
			vis[ i * prime[ j ] ] = 1;
			if( i % prime[ j ] == 0 ) {
				phi[ i * prime[ j ] ] = phi[ i ] * prime[ j ];
				break;
			}
			else
				phi[ i * prime[ j ] ] = phi[ i ] * ( prime[ j ] - 1 );	
		}
		f[ i ] = ( f[ i - 1 ] + 1ll * i * i % p * phi[ i ] % p ) % p;
	}
}

map< long long , int > Map;
int sum2( long long x ) {
	x %= p;
	return 1ll * x * ( x + 1 ) % p * ( 2 * x + 1 ) % p * inv6 % p;
}
int sum3( long long x ) {
	x %= p;
	return 1ll * x * x % p * ( x + 1 ) % p * ( x + 1 ) % p * inv4 % p; 
} 
int Sumphii( long long n ) {
	if( n <= MAXN ) return f[ n ];
	if( Map[ n ] ) return Map[ n ];
	
	int Ans = sum3( n );
	for( long long l = 2 , r ; l <= n ; l = r + 1 ) {
		r = n / ( n / l );
		Ans = ( Ans - 1ll * ( sum2( r ) - sum2( l - 1 ) ) % p * Sumphii( n / l ) % p ) % p; 
	}
	return Map[ n ] = ( Ans + p ) % p;
}
int solve( long long n ) {
	int Ans = 0;
	for( long long l = 1 , r ; l <= n ; l = r + 1 ) {
		r = n / ( n / l );
		Ans = ( Ans + 1ll * ( Sumphii( r ) - Sumphii( l - 1 ) ) * sum3( n / l ) % p ) % p;
	}
	return ( Ans + p ) % p;
}

int Quick_pow( int x , int po ) {
	int Ans = 1;
	while( po ) {
		if( po & 1 ) Ans = 1ll * Ans * x % p;
		x = 1ll * x * x % p;
		po /= 2;
	}
	return Ans;
}
int Inv( int x ) {
	return Quick_pow( x , p - 2 );
}

signed main( ) {
	scanf("%d %lld",&p,&n);
	sieve( ); 
	inv4 = Inv( 4 ) , inv6 = Inv( 6 );
	printf("%d", solve( n ) );
	return 0;
}